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INTRODUCTION 


Nonconformal contact machine elements in power train systems such as 
gears, rolling element bearings, and cam and follower mechanisms are subject 
to transient lubrication. The transient characteristics are due to the time 
variation of loading, geometry, and the rolling or sliding speed in the line 
or point contact. These variations result in a squeeze effect which affects 
the minimum film thickness distribution. An example of this is the ball 
bearings in a rotordynamic system in which there exist cyclic variations of 
the dynamic load. Recently, the transient hydrodynamic and elastohydrodynamic 
line contact problem has received much attention (Refs. 1 to 3). Among the 
several authors, Vichard (Ref. 1) pioneered the basic transient 
characteristics of the line contact problem analytically and experimentally 
Including the viscous damping phenomenon. In this paper, the transient 
solution of the hydrodynamical ly lubricated point contact presented. 

In solving the point contact transient problem numerically, a fast 
computer code is needed to solve the two dimensional Reynolds equation for 
many time steps. Numerical methods for solving the simultaneous equations 
resulting from the discretization of the Reynolds equation are usually 
performed using either iterative methods or semidlrect methods (Ref. 4). The 
former commonly involves the Gauss-Seidel method, the latter combines the 
Newton-Raphson method with a direct inversion of the Jacobian matrix. An 
Important difference between the iterative method and the semidlrect method is 
that the initial guess plays an important role in the latter, whereas the 
former Is relatively insensitive to the initial guess. With the semidlrect 
method, the use of a previous solution as an initial guess accelerates the 
solution process, but a good initial guess usually does not help the iterative 
method significantly (Ref. 4). The semidlrect method is preferred for 
transient problem since the solution of the previous time step accelerates the 
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ABSTRACT 

The transient analysis of hydrodynamic lubrication of a point-contact is 
presented. A body-fitted coordinate system Is Introduced to transform the 
physical domain to a rectangular computational domain, enabling the use of the 
Newton-Raphson method for determining pressures and locating the cavitation 
boundary, where the Reynolds boundary condition is specified. In order to 
obtain the transient solution, an explicit Euler method is used to effect a 
time march. The transient dynamic load is a sinusoidal function of time with 
frequency, fractional loading, and mean load as parameters. 

Results Include the variation of the minimum film thickness and phase-lag 
with time as functions of excitation frequency. The results are compared with 
the analytic solution to the transient step bearing problem with the same 
dynamic loading function. The similarities of the results suggest an 
approximate model of the point contact minimum film thickness solution. 


*NASA Resident Research Associate at Lewis Research Center. 



next step solution. Furthermore, the Newton-Raphson method has a quadratic 
convergence rate, so, In general, the solution can be terminated within ten 
Iterations. When a parallel processing computer using vectorlzatlon Is 
employed the matrix Inversion is very fast. In addition, there Is no need to 
use underrelaxatlon factors, and the solution can be obtained more rigorously 
than Is typical with Iterative methods. The matrix Inversion can be done by 
the Thomas algorithm, and there Is no need to store the whole Jacobian matrix. 

When the semldlrect method Is used In the point contact problem, the 
cavitation boundary, where the Reynolds boundary condition (B.C.) Is specified, 
Is difficult to locate. There Is a fundamental difference between the line 
contact and the point contact problem. In the line contact case, the Reynolds 
equation Is Integrated once; the Neumann condition is Introduced; and the 
Integration constant Is found as a part of the solution. In the 
two-dimensional problem, the Reynolds equation can not be Integrated. Since 
the Reynolds B.C. Insures mass conservation across the boundary, the 
cavitation boundary should be located as accurately as possible. However, the 
location Is not known In advance; it is a part of the solution. It Is a free 
boundary where two B.C.'s are present: Dlrichlet B.C. (pressure Is zero), and 
Neumann B.C. (normal pressure gradient Is zero). The relaxation method of 
Chrl stopherson (Ref. 5), derived for the hydrodynamic lubrication of a journal 
bearing, has been used to solve this kind of free boundary value problem. 

This method truncates negative computed pressures whenever they occur during 
Iteration. However, this method can not be used In the semldlrect method. In 
this work a body-fitted coordinate system Is Introduced which transforms the 
unknown boundary Into a fixed boundary and the unknown boundary function Is 
Introduced Into the equations of motion. The smooth cavitation boundary is 
found up to truncation and machine errors, whereas the result for 
Chrl stopherson ' s method Is dependent upon the mesh size near the boundary. To 
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detect the minute change of the cavitation boundary between the adjacent time 
steps, the current method Is desirable. Another advantage of this method Is 
that a nonzero pressure gradient condition can be implemented for very lightly 
loaded cases where surface tension may play an important role, or for 
non-Newtonian, viscoelastic fluids. 

In the present paper the transient hydrodynamic lubrication of a step 
bearing Is solved analytically to provide physical Insight into the transient 
characteristics of hydrodynamic lubrication. Next, the point contact problem 
Is solved numerically by the Newton-Raph son method with Thomas algorithm. 

This method is fast and does not require vast computer storage. Parallel 
processing by vectorlzatlon is also utilized. 

The variation with time over a loading cycle of the minimum film 
thickness, squeeze velocity, and the cavitation boundary Is studied for a wide 
range of excitation frequencies. 

NOMENCLATURE 

F dimensionless load 

Fq dimensionless mean load 

F right hand side equation of discretized equation 

f load, N (point contact), N/m (step bearing) 

fQ mean load, N (point contact), N/m (step bearing) 

G dimensionless cavitation boundary function 

G 1 first derivative of G with respect to Y 

G" second derivative of G with respect to Y 

g cavitation boundary curve function 

H dimensionless film thickness 

Hq dimensionless minimum film thickness 

Rg normalized dimensionless minimum film thickness, Hq/Hq,,, 
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Hom dimensionless minimum film thickness for Fq 
h fl 1m thickness, m 

ho minimum film thickness, m 

k number of Iteration of Newton-Raphson method 
L length of the step bearing, m 

% reference length for order-of-magnl tude analysis, m 
NI number of grid In \ direction 
NJ number of grid In n direction 

n normal direction vector 

P dimensionless pressure 

p pressure, N/m 2 

R radius of sphere, m 

R residual vector of discretized equation 

t time, sec 

u solution vector of the discretized equations 
u m average surface velocity In x-dlrectlon, m/sec 
uq reference velocity for order-of-magnltude-analysls, m/sec 

X dimensionless coordinate along rolling direction 

X/\ dimensionless Inlet boundary location In X-dlrectlon 

x coordinate along rolling direction 

x/\ Inlet boundary location In x-dlrectlon 

Y dimensionless coordinate transverse to rolling direction 
Y3 dimensionless Inlet boundary location In Y-dlrectlon 
y coordinate transverse to rolling direction 

yg Inlet boundary location In y-dlrectlon 

a viscosity-pressure coefficient, m 2 /N 

a dimensionless viscosity-pressure coefficient 
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0 fractional loading amplitude for sinusoidal loading 
Y dimensionless frequency 

6 dimensionless film thickness of the step bearing 

5 normalized film thickness of step bearing, S/SQm 

$0m dimensionless film thickness of the step bearing for mean load 

n lubricant viscosity. Pa-sec 

p dimensionless lubricant viscosity 

po lubricant viscosity at atmospheric pressure, Pa-sec 

v kinematic viscosity, m 2 /sec 

C,n coordinates of transformed domain 

p lubricant density, kg/m^ 

t dimensionless time 

<t> s phase angle of the step bearing solution, deg 

4>p phase angle of the point contact solution, deg 

Q physical domain 

O' computational domain 

<o frequency of sinusoidal loading, (cycle)/sec 

ANALYTICAL SOLUTION OF A STEP BEARING 
Consider the simple step bearing shown In Fig. 1. Note that the step 
bearing used here is subjected to an oscillating normal motion and is closed 
at the exit end. To the authors' knowledge, this particular solution Is not 
available In the literature and is therefore presented here. The film profile 
and the dynamic force are: 

h(x , t) = h(t) , 0 < X < L, 

- 0, x « L, (1) 

f(t) = fo<l + 0 sin wt). (2) 
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For an Incompressible, Isovlscous, Newtonian fluid, the governing 
equation Is, 


lx ( h3 

3fi\ 
a xy 

■ ,2 V b 

9h 

8x + 

,, 9h 
% 8x 

U 1 

V -2* 

(3) 

The boundary conditions 

and 

the 

Initial condition 

are, 




P 

= 0 

at 

x = 0, 

J 




h 

= 0 

at 

x = L, 


(4) 
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when 
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the dimensionless equatl 

ons 
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. X <1 , 





= 

o, 

X = 1, 


(5) 


P = 0 at X = 0, 

6 = when x = 0, 


F(t) ■ Fq(1 + p sin yx) • (6) 

After Integrating Eq. (5) three times using, 

f 1 P dX = F(x> (7) 

J 0 

the following nonlinear first order differential equation results: 


2_ 8S_ 3_ 

S 3 8t “ s 2 


1 

2 


F(x>. 


The solution of Eq. (8) subject to the B.C.'s in Eq. (5) Is, 


( 8 ) 


S(x) 


J_ _ ^0 F 0 A -3t !o V 
s 2 6 2<y 2 + 9)/ 6 + 2(y 2 + 9) 


,-1/2 


(3 Sin yx - y COS yx) 


< 9 ) 

After a sufficiently long time (x -> ®), the exponential term vanishes to zero, 
and the time variation of the film thickness becomes. 
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1 - 1/2 


S(x) = 6 


Om 


where 


Om 


V"^ 

Vt- 


s1n(yx - 4* s > 


♦s = tarr1 ©* 


( 10 ) 


The formula for the squeeze velocity Is obtained by differentiating 
Eq. (10), 


om 


- (- 8 


JL 


'{ 7 . 


COS(yx - <J> ) 


-,-3/2 


"{ 7 . 


sin(yx - 4> s ) 


( 11 ) 


ANALYTICAL FORMULATION OF THE POINT CONTACT PROBLEM 


The physical model is Illustrated in Fig. 2. The radius of the sphere is 
R and the dynamic force Is the same as that of the step bearing. The two 
dimensional, transient, Incompressible form of the Reynolds equation for 
Newtonian flow Is, 




where 


p = p(x,y ,t) 
h * h(x,y,t) 


( 12 ) 


V - p(x,y,t>. 

The parabolic approximation of the film thickness equation of the sphere 


Is: 


h = h 0 + IR (x2+ y2) 


(13) 


At a given time, the generated pressure distribution Is balanced by the 
dynamic load, 


f(t) 


P(x,y,t)dx dy. 


(14) 
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The plezovl scous effect Is modelled by the Barns relation (16): 


« a P 

v = V • 

The boundary conditions are: 

P - 0 at x = x A 0 < y < y g , 
p = 0 at x A < x < g(y B ,t) y = y B , 

p = 0; ^ = 0 at x = g(y,t) 0 < y < y R 

3n 0 

* 0 at X A < x < g(0,t) y = 0. 

At the cavitation boundary, x = g(y,t), the pressure and the normal 
pressure gradient are zero (Reynolds B.C.). Using symmetry at y = 0, the 
Neumann B.C. Is Imposed and only half of the domain Is modelled. 

With the following definitions, 



X - 


x . 
R’ 




P _ _&B_. 


F = 


f 

•W 


x 



Y 


the dimensionless equations are: 





f h£ 3P ) ^ a_ f h£ 3P ] 3H ,, 3H 
3X Vj} 3X7 + 3Y Vjl 3y J = 12 3X + 12 fr’ 

H - H Q + i (X 2 + Y 2 ), 


F(x> = 2 


„Y b p G(Y) 


0 J 


P(X,Y,x)dX dY, 


F(t) = F n (l + p sin y t) , 


P = e 


aP 


-8 


with a * 1.5131X10 In this study. 


( 15 ) 


(16) 


(17) 
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To fix the unknown cavitation boundary, the following body-fitted 


coordinate transformation shown in Fig. 3 is Introduced: 


£ 


v x - V 

G(Y,t) - X A ’ 


n = Y, (18) 

| J | = (G(Y,t) - X A )/Yg. 

|J| is the Jacobian of the coordinate transformation which shows that as 
long as G(Y.-t) Is not equal to Xy\, there exists a conformal mapping between 
the physical domain and the computational domain. 

The differentiations transforms to the following: 

Y 

3_ B 8_ 

9X = G - X A 3C’ 

3_ d_ EG' _3_ 

3Y = 8 n " G - X A 3£’ 


2 2 2 ’ 
3X^ (G - X A )^ dC 


2 2 2 2 
_3_ r(G’r 3 2CG * 


*[« 


3Y 


2 (G - X.) 2 3E 2 G - X A ^ 


8- 3‘ -L 2<G ') 2 - G " <S -V] 

+ n + ^ 


, A , 3r A 3n 

The Reynolds equation In the (£,n) system Is, 


(G - X A )‘ 


where 


AP CC + BP C + CP + DP_ + EP + F = 0, 
££ £n nn £ n 


A . A 3 [y 2 * 5 2 <G') 2 ], 
B = A 3 C-2?G'(G-X a )], 

C = A 3 (G - x A ) 2 , 


a£‘ 


(19) 


(20) 


D = A 1 Yg(G - X A ) - A 2 5G'(G - X A ) + A 3 c[2(G') 2 - G"(G - X ft )] , 
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E = A 2 (G - X A )\ 


F = -A 4 (G - X A ) , 


. 3H_ 9H „3 3_ / 1 \ 

1 = - 3X + n 3X - ’ 


a _ IMI M H 3 3_ /I] 

2 “ - 3Y + " 3Y 1 - r 


A - — 


A 19 9H 19 9H 

a 4 ■ 12 ax * 12 H' 

In the above formulation, A], A£, A 3 , and A 4 can be transformed to the 
(£,n> coordinate system using Eq. (19). At the cavitation boundary. 


3P 

3n 


1 


* J77 


,2 L 


(G 1 ) 

Since 3P/3 n = 0 at 5 = Y g 


1 /v c , ri ,2.3P r.3P 

G - X A <Y B + ?<G } > 3^ " G 3 n 


= 0. 


( 21 ) 


9P „ 0 

3E ’ 


P - 0 at £ » Yg. 


At n = 0 (Y = 0), the symmetry condition Is, 

3P _ 3P _ EG' 3P n 
3Y " 3 n G - X A 3E " 


( 22 ) 


(23) 


But, G 1 = 0 due to the symmetry of cavitation boundary and It follows 

that, 

P = 0 at n = 0. (24) 

The transformed film thickness equation and the force balance equation 
are expressed, 
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H<S,n> = H q + 


F(t) 



5<G - 


X A> 


-i 2 


+ n 


■ y B (G - x a } 

B P<S,n,T> — y-5- dl- d n . 

0 t B 


(25) 


(26) 


In the above formulation, the unknown boundary curve function G Is 
introduced Into the governing equations while the computational domain Is 
fixed. 

NUMERICAL METHODS 

Equation (20) Is a nonlinear partial differential equation. The 
nonlinearity is due to the piezoviscous relation and to the function G In 
the transformed Reynolds equation. 

Spatial Discretization 

In order to minimize the number of grid points while maintaining accuracy, 
a smoothly varying nonuniform spacing is generated by a two-sided stretching 
function, (hyperbolic tangent) (Ref. 10). The finest spacing is near the 
cavitation boundary which Is also near the maximum pressure gradient. 

Figure 4 shows the finite difference mesh structure. The Increments in 
S and n and are such that 


*1 - * 1-1 


AC 


*1+1 ‘ * 


j * 


(27) 


n-i - n 


J-l 


= An 


n J+l ” n J " r n Ari ' 

By the Taylor series expansion, the finite difference approximations of 
derivatives with respect to § and n are, 
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3P _ r g P 1-1, J + (r E ' 1)P I,J + P I+1.J 


as 


r^(l + r^)A£ 


p -r 2 P T + (r 2 - 1)P T . + P T . . 
oP _ n I , J-l n I , J I , J+1 

3n r (1 + r )An 

n n 

a^P = 2 r S P I-l,J ~ (r S + 1)P I.J + P I+1 . J 
a ? 2 r^(l + r ^)^ 2 

£2 , 2 r n P I,J-1 ~ (r n * 1)P I.J + P I.J+1 


3n 


2 

3^P 


r (1 + r )An 
n n 


3£3n r c r (1 + r c )(l + r )A£An 
s n $ n 


r e r n P I-l , J-l ■ r n <r 5 ■ 1 >P I , J-l ■ r n P I+l .J-l ‘ r 5 <r n ' ,)P I-1,J * <r f ‘ n 


X <r n - ' )P I.J * <r n - ,)P W.J - r £ P I-l ,J+1 + <r C - ,)P I.J*. * P I+1 ,J+1 


(28) 

Substituting Eq. (28) Into the transformed Reynolds Eq. (20) the following 
discretized equations results, 

R I,J = C 1 P I-1,J+1 + C 2 P I , J+l + C 3 P I+1 , J+l + C 4 P I— 1 , J + C 5 P I , J + C 6 P I+1 , J 

+C 7 +CP +CP + C = 0 , (29) 

7 I-1,J-1 8 I, J-l 9 I+l.J-l 10 


wl th 


P 1 ,J = P NI , J " 0 1 < J < NJ, 


P I,NJ " 0 
P I,0 = P I ,2 
P NI + 1 ,J = P NI-1 ,J 


1 < I < NI, 

1 < I < NI, 

1 < J < NJ. 


13 



Steady-State Solution Method 

The transient solution Is formed by computing the steady-state solution 
for each time step Including the squeeze term. The numerical technique for 
the steady-state solution along with the Thomas algorithm and Newton-Raphson 
method Is described first. 

The discretized form of transformed equation Is, 

K(u)u = F (30) 

The vector u represents the unknown values, pressures and cavitation 
boundary. For an Isovlscous condition K(u), contains the function G, and, 
for a plezovlscous condition, It Includes pressures as well. The discretized 
simultaneous equations are nonlinear. Even for the linear free boundary value 
problem. It has a nonlinear characteristics since the unknown boundary Is 
associated with the solution. 

The Newton-Raphson method Is described, 

vi ■ 5 k - J ''‘V S< V <3,) 

where R(u^) = K(u k )u k " ^ res ^ual vector and J < u ) Is the Jacobian 

of the system of equations. In practice, the Iteration Is organized as, 

, .4 , 4 f 4 . -4 -4 -4 

JCu^Au^ = -R(u k ), u k+1 = u k + Au^. (32) 


For this study, the vector u Is, 

3 - ( P z,j- P 3.J P NI-1 .J- G j) T - 3 <33) 

In which P] j and Pn^j are zero from the Dlrlchlet boundary condition. 

-4 

The residual vector R Is, 


R 



R NI-1 ,J 



J = 1, NJ - 1 


(34) 
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The Jacobian matrix is a block tridiagonai matrix In Fig. 5, and each 
block Is a one-sided arrow-shaped matrix, Fig. 6. In the formulation of each 
block matrix of the Jacobian, the last columns are the differentiations of the 
residual vector with respect to the cavitation boundary function, G. Since all 
the coefficients in the discretized Reynolds equation are composed of Gj, G'j, 
and G ’ 1 j , it is easier to calculate them numerically (Ref. 11) using: 


9R 


I.J 1 


9G- 


R I,J (G J + c 


g ,w i,j ) 


r i,j (g j 


w i,j ) 


(35) 


where W I,J contains all other variables except Gj. The value of 
eg can be chosen to be sufficiently small not only to maintain good accuracy 
of Eq. (35) but to prevent serious round-off errors. In this calculation, 
eg is set to 10 - 9 in double precision. 

The block tridiagonal system of Eq. (30) is solved by the Thomas 
algorithm (Ref. 12). This algorithm inverts the whole matrix at a time by 
matrix multiplication and inversion of the block matrix, which is quite fast 
on a parallel processing computer with small memory storage size equal to 
2 x NI x NJ x NJ. The matrix inversion is accomplished using UNPACK. 

The Newton-Raphson method requires a good Initial guess of the solution. 
For this purpose, the Gauss-Seidel iteration method is used to get an 
approximate pressure distribution and cavitation boundary location. Once one 
solution is obtained by the Newton-Raphson method, it is used for the guess to 
next solution. The convergence criteria are 


(1) pressure 



< 1.0x10 


-4 
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(2) cavitation boundary 



1 .OxlO -4 


(3) force balance 

^ input " ^output . 

= < 1.0x10 

Input 

In order to make sure of the convergence, the L 2 -norm of the residual 
vector Is monitored. In general, the solution converges within 3 to 8 
iterations. In this study, NI - 41 , NJ = 31 . 

Transient Solution Method 

For the steady-state solution, the problem is to find Hq for a given 
load, or for a hydrodynamic case, the load capacity can be calculated for a 
given Hq. But, for the transient case, there Is an additional unknown value 
to be determined, the squeeze velocity. The basic solution technique Is to 
use a "time-march." That is, Hq is fixed from the previous time step, and 
the squeeze velocity is found that balances the generated pressure distribution 
with the dynamic force at that time. The detailed computation procedure is 
provided In Fig. 7. At the first time step, the steady-state Reynolds 
equation is solved to find Ho m . and, fixing Hq, the transient Reynolds 
equation is solved including the squeeze term to find the squeeze velocity 
using the force balance equation. For this purpose, a bisection method Is 
used, with an approximate range of squeeze velocities according to the history 
of dynamic force and the minimum film thickness variation. Once a converged 
solution is obtained, the minimum film thickness of the next time step Is 
estimated from the following expression: 
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(36) 


C 1 = Hq + Ax, (n = present time step) 

The film thickness and squeeze velocities are established at successive 
time steps and the calculation Is continued beyond the first complete loading 
cycle until the periodic requirement Is reached. The convergence criterion Is: 

_ 4 

< 1.0x10 , IC = number of cycle 

In this calculation, 361 time steps with 1° Increment are used in one 
loading cycle. 

RESULTS AND DISCUSSION 

The analytical solution of the step bearing demonstrates that 6 approaches 
one with a phase-lag of 90° as y increases (Fig. 8). This asymptotic 
behavior Is due to the squeeze action caused by the dynamic forces. Figure 9 
shows the squeeze variation of Eq. (11). This phenomenon Is physically 
similar to a nonlinear massless spring-damper system with forced vibration 
shown in Fig. 10, sometimes referred as a "half a degree of freedom system." 

The response of this system is that the amplitude approaches a constant value 
and the phase-lag goes to 90°. Although the transient solution of the point 
contact problem can not be solved analytically and requires numerical 
computation, it may be speculated that basically It also has a similar 
nonlinear spring-damper system. In the following example, the numerical 
results of the point contact problem are compared to the step bearing solution. 

For this study, Fq = 3000 and (3 = 0.3 with different y's. The 
minimum film thickness for Fq is 1 . 247 lxl 0 — ^ for the isoviscous case and 
1 . 3907x1 0 - 5 for the piezoviscous effect with X/\ = 0.08, Yg = 0.06. Figure 11 
shows the pressure distribution for Fq and Fig. 12 delineates the detailed 


K’n 0 *' - (H 0>n C 

<vJ c 
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cavitation boundary curve in which the minimum value of G occurs at Y = 0 
and it Increases up to a certain location and then decreases because of the 
geometry of the sphere. 

Figure 13 illustrates the time variation of the normalized minimum film 
thickness (Rq> during one loading cycle with 361 time steps. The squeeze 
velocity distribution is shown in Fig. 14 for different y's. These results 
are qualitatively similar to those of the step bearing solution. However, it 
should be noted that the order of the nondimensional excitation frequencies is 
different since L is used as a reference length in the step bearing while 
R is used for the point contact case. 

Equation (10) may be put in the following form. 


6 = 


11/2 


1 

1 + a $ |3 si n(yx - <t> $ ) 


where 



1 

1 + (X $Y ) 2 


4> s = tan -1 (X $ y) 



(37) 


The variation of a s and 4> s are plotted in Fig. 15. 

For quantitative analysis of the transient point contact problem, the 
following formula is suggested by Eq. (37), 


H 0 


.1 + a p |3 sin(yT - <|>p 


»] 


(38) 


Equation (38) is deduced based on the fact that Hom is Inversely 
proportional to Fq whereas $Q m to ^Fq- The unknown values in Eq. (38), 

ap and 4>p, are obtained by a nonlinear least square fit with 361 data 
point. Figure 16 shows the comparison between the numerical results and the 
curve fit. The best curve fit can be obtained by letting the numerator of 
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Eq. (38) be variable, however, It Is near 1, for example, 1.005 for y = 100, 
1.019 for y = 1000. The curve fitting results are recorded In Table 1. 

Fig. 17 shows the variation of a p and $ p , qualitatively, similar 
characteristics to the analytical step bearing solution with different order 
of magnitude of y (Fig. 15). The value of X p Is obtained assuming the 
following relation, 

<j> p = tan -1 (X pY ) (39) 

X p is nearly constant over a wide range of y, approximately 0.0054. If an 
analytical solution were possible, the a p would be a function of X p . 

However, since it also would be a function of the geometry associated with the 
cavitation boundary, no attempt is made to obtain a form similar to Eq. (37). 
Instead, for design purposes, Eq. (38) can be used along with Table 1. 

For the plezoviscous solution, a p is smaller than that of the isovlscous 
solution (Fig. 17), but <j> p ' s are virtually the same. The a p 's 
asymptotically approach those of the isoviscous case. Figure 18 shows this 
more vividly. Due to the plezoviscous effect, the distribution of Rq is 
more damped with the same phase angle. The X p 's for the plezoviscous case 
are nearly constant and equal to the isoviscous case (see Table 1). This 
implies that X p is a characteristic of the transient point contact problem 
of the current model. 

Figure 19 illustrates the location of the outlet boundary at Y = 0 
normalized by that for the steady-state solution of Fn. For the steady state 
case, G(0,x) approaches the point of contact as the load increases. However, 
when y is greater than zero there exists a substantial variation in G(0 ,t) 
due to the squeeze action. When the squeeze is downward, G(0,x) may be 
stretched outward and vice versa. For example, when y = 200, there is a 
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downward action between a-b and c-d in Fig. 19, and upward action between 
b and c. These points correspond to those In Fig. 14. 

In the foregoing analysis, the Reynolds Eq. (12) neglects the inertia 
forces. But, as y Increases, the validity of this assumption becomes 
suspect. This assumption Is examined by an order-of-magnltude-analysl s 
of the steady-state Navler-Stokes equation In Ref. (12). When the modified 
Reynolds number Is much less than one. 



(40) 


the Inertia forces can be neglected. Here, ug Is a reference velocity, 
2. Is a reference length In the x-dlrection, and hg Is that in the film 
thickness direction. Using, 


U 0 = Rco; y = jp~ <41> 

m 

the following relation for the validity of the assumption that Inertia 
forces are negligible Is, 


For example, If Hg = 10 - 5, R 
v = 10 - 5 m 2 /sec, 


h A h n u m 

H - — Rp - — m 

H 0 ' R ’ Re " v 1 


10 -2 m, u = 0. 1 m/sec, and 


(42) 


Y « 10 8 (43) 

Even for y = 1000, Inertia effects remain negligible. 


CONCLUSIONS 

The transient solution of the hydrodynami cal ly-lubrlcated point contact 
problem Including the squeeze effect Is obtained numerically using the 
ball-on-plane model, A new computational algorithm Is Implemented to deal 
with the cavitation boundary by the semidirect method with the advantage of 
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supercomputi ng. This method provides a faster and more rigorous way to solve 
the nonconformal contact problem with a Newtonian fluid than the conventional 
iterative method, and the flexibility to deal with more complex boundary 
conditions for lightly loaded bearings and more realistic rheological models. 

The qualitative and quantitative analysis is compared with the analytical 
solution of a dynamically loaded step bearing solution using a nonlinear curve 
fitting method. It is found that there exists a characteristic similarity In 
the transient responses to a nonlinear massless (i.e., no Inertia) 
spring-damper system, in terms of the variation of the minimum film thickness 
and phase angle. According to an order of magnitude analysis, it is confirmed 
that the inertia-forces are negligible for a wide range of practical 
excitation frequencies. 

These results can be applied to the design of moderately loaded ball 
bearings In rotordynamic systems and can be extended to gear design adding the 
time variation of the geometry and speed. For highy loaded elliptical contact 
case, the elastic deformations and ellipticlty parameter need to be considered. 
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TABLE 1 - CURVE FITTING RESULTS OF EQ. (38) 


Isovi scous 


Y 


D 

X P 

■a 

D 

X P 

0 

0.890 

0.0 


0.734 

0.0 


25 

.882 

7.0 

0.00491 

.728 

7.0 

0.00491 

50 

.859 

14.7 

.00523 

.709 

14.7 

.00526 

100 

.781 

28.3 

.00538 

.646 

28.4 

.00541 

150 

.687 

39.2 

.00544 

.571 

39.2 

.00544 

200 

.599 

47.5 

.00545 

.499 

47.4 

.00544 

250 

.523 

53.8 

.00546 

.438 

53.7 

.00545 

300 

.461 

58.7 

.00547 

.386 

58.5 

.00544 

350 

.409 

62.3 

.00545 

.343 

62.3 

.00545 

400 

.368 

65.3 

.00544 

.309 

65.2 

.00541 

500 

.303 

69.8 

.00544 

.255 

69.8 

.00544 

750 

.210 

76.1 

.00539 

.176 

76.1 

.00539 

1000 

.159 

79.4 

.00536 


79.4 

.00536 


/{ t) =/o(1 “Bsinojt) 



Figure 1. - Schematic view of the step bearing configuration. 
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Figure 2. ♦ Physical model of the point contact problem. 
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(a) Physical domain. (b) Computadonai domain. 

Figure 3. - Coordinate transformation of the physical domain to the 
computational domain. 



Figure 4. - Finite difference mesh structure. Figure 5. - The Jacobian matrix of eg. [32]. 
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Figure 6. - The elements in diagonal block matrix. 



Figure 7. * Flow digram for the transient solution. 
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TIME STEP 

Figure 16. - Comparison of the numerical results with 
curve fit for 100. 
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